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Abstract. Precession driven flows are found in any rotating container filled with 
liquid, when the rotation axis itself rotates about a secondary axis that is fixed in 
an inertial frame of reference. Because of its relevance for planetary fluid layers, 
many works consider spheroidal containers, where the uniform vorticity component 
of the bulk flow is reliably given by the well-known equations obtained by Busse in 
1968. So far however, no analytical result on the solutions is available. Moreover, the 
cases where multiple flows can coexist have not been investigated in details since their 
discovery by Noir et al. (2003). In this work, we aim at deriving analytical results on 
the solutions, aiming in particular at, first estimating the ranges of parameters where 
multiple solutions exist, and second studying quantitatively their stability. Using the 
models recently proposed by Noir & Cebron (2013), which are more generic in the 
inviscid limit than the equations of Busse, we analytically describe these solutions, their 
conditions of existence, and their stability in a systematic manner. We then successfully 
compare these analytical results with the theory of Busse (1968). Dynamical model 
equations are finally proposed to investigate the stability of the solutions, which allows 
to describe the bifurcation of the unstable flow solution. We also report for the first 
time the possibility that time-dependent multiple flows can coexist in precessing triaxial 
ellipsoids. Numerical integrations of the algebraic and differential equations have been 
efficiently performed with the dedicated script FLIPPER (supplementary material). 
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A rotating solid object is said to precess when its rotation axis itself rotates about a 
secondary axis that is hxed in an inertial frame of reference. In this work, we consider 
a precessing spheroid of fluid, in rotation around its symmetry axis. This geometry is 
indeed relevant for planetary fluid layers, such as the Earth liquid core ca, where 
precession-driven flows may participate in the dynamo mechanism generating their 
magnetic helds HI El [23]. These flows may also have an astrophysical relevance, for 
instance in neutron stars interiors where they can play a role in the observed precession 
of radio pulsars [9]. Finally, precessing spheroids have been studied as turbulence 
generators, in particular when the angle between the container rotation axis and the 
precession axis is 90° mm- 

The hrst theoretical studies of this spheroidal geometry considered an inviscid fluid 
[m [30l 1^ . Assuming a uniform vorticity, they obtained a solution, called Poincare 
flow, given by the sum of a solid body rotation and a potential flow. However, the 
Poincare solution is modihed by the apparition of boundary layers, and some strong 
internal shear layers are also created in the bulk of the flow m^- In 1968, Busse have 
taken into account these viscous effects as a correction to the inviscid flow in a spheroid, 
by considering carefully the Ekman layer and its critical regions (see also lasiilla]). 
Based on these works, [6] and iza have proposed models for the flow forced in precessing 
triaxial ellipsoids. Note that, beyond this correction approach, the complete viscous 
solution (including the hne description of all the flow viscous layers) has recently been 
obtained in the particular case of a weakly precessing spherical container HT]. 

When the precession forcing is large enough compared to viscous effects, instabilities 
can occur, destabilizing the Poincare flow. First, the Ekman layers can be destabilized 
[20] through standard Ekman layer instabilities [T^ [8]. In this case, the instability 
remains localized near the boundaries. Second, the whole Poincare flow can be 
destabilized, leading to a volume turbulence: this is the precessional instability [23] . 
This small-scale intermittent flow conhrm the possible relevance of precession for energy 
dissipation or magnetic held generation, and has thus motivated many studies. Early 
experimental attempts [lOl [39] to conhrm the theory of Busse [5] did not give very 
good results [28] . Simulations have thus been performed in spherical containers |32l EE] , 
spheres [26], and hnally in spheroidal containers [211122], allowing a validation of the 
theory of Busse [5]. Experimental conhrmation of the theory has then been obtained 
in spheroids [ 23 ], a work followed by many experimental studies involving spheres 
[iniiiHiE], and spherical containers [38] . 

Finally, the dynamo capability of precession driven hows has been demonstrated in 
spheres [331E1], spheroids [31] and cylinders [27], allowing the possibility of a precession 
driven dynamo in the liquid core of the Earth [16] or the Moon [7]. 

In this work, we focus on the precession forced how described by the system of 
algebraic equations obtained by Busse in 1968 [5] for precessing spheroids. As shown 
by Noir and co-workers [23], these equations reliably describe the how but can lead to 
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Figure 1. Sketch of the problem under consideration. We consider a container filled 
with liquid and set in rotation at tio along its symmetry axis Oz. The spheroid axis 
is tilted at the precession angle a and fixed on a rotating table, which rotates at the 
precession rate ftp. In a planetary context, the fluid corresponds to a liquid core, the 
container to the mantle, and the rotating table plane to the orbital plane, i.e. the 
ecliptic for the Earth. 


multiple solutions for particular ranges of parameters . So far, these multiple solutions 
cases have not been investigated in details. Indeed, the study of Noir and co-workers [21] 
is the only work considering these possible multiple solutions, and they only calculate 
the stability of the solutions in certain cases. In particular, the ranges of parameters 
allowing multiple solutions are not known (which is also partly due to the absence of 
any analytical result). No analysis of these equations have been performed to obtain 
rigorous constraints on the solutions, or to estimate the solutions and their stability. 

We propose here to tackle analytically these issues in order to obtain analytical 
estimates and scaling laws, to compare our results to the exact solutions, and to 
investigate analytically the solutions stability. 

In section [21 we introduce the problem considered in this work in a general 
framework. In section [3l we hrst present few multiple solutions cases (section 13.11) 
and then, in section 13.21 we introduce recently proposed theoretical models [25], 
which extends the Busse equations into a dynamical framework. Relying on this 
new approach, theoretical investigations are tractable (section and the obtained 
analytical estimates are then compared with the predictions of the Busse model. The 
stability of the solutions are studied in section 0] 

2. Mathematical description of the problem 

We consider an incompressible homogeneous fluid of density p and kinematic viscosity 
u enclosed in a spheroidal cavity, rotating along its symmetry axis at fio = (^z 

is a unit vector), and processing at the angular velocity flp, with a precession angle 
a , also called obliquity (see hgureH]). Noting c the length of the spheroid symmetry 
axis, we use the length a of the other principal axis as the length scale. Using l/VLp 
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as a time scale, the problem is completely defined by four parameters: the oblateness 
rjs = 1 — c/a, or equivalently the flattening ellipticity e = (c^ — a^)/(c^ + a^), the Ekman 
number E = the precession angle a, and the Poincare number P = VLp/VLo 

in the literature. The dimensionless components of the precession vector are then 
naturally Px = Psina, Py = 0, and P^ = Pcosa (the axisymmetry naturally allows 
the choice P^ = 0 without any loss of generality). For planetary applications and 
turbulence generators, works of the literature typically consider moderate values of e, 
small Ekman number P -C 1, and values of a are typically moderate for planetary 
cores {a ~ 23.4° 0.41 rad for the Earth), or large for turbulence generators studies 

(a = 7r/2) [iniin]. 

Considering the dimensionless fluid rotation rate fio, which is half the dimensionless 
vorticity, one can derive a system of equations governing the uniform vorticity bulk 
component of the flow (e.g. mn) 

nl + nl + nl-n, = o, (i) 

-P, Vty = r /3 VtyVt, + {\r ^x^y^ + Ai nyn-^/y^^E, (2) 

P.^y = - \r (1 - f^.) ^/P, (3) 


which are exactly equations (20)-(22) of [21], or the equations (21)-(23) of |6] in the 
particular case of a spheroid (772 = 0 in their notations, and Py = 0). As shown by 
|24j . this system of equations is equivalent to the well-known implicit expression (3.19) 
of Busse [ 5 ]. Equation o is the so-called no spin-up condition (see the solvability 
condition 3.14 of [5], or equation 12 of [21]) given that it imposes (f2 — efl) ■ = 0, i.e. 

it forbids any differential rotation along 17. Equations ([2])-(l3]) are simply obtained from 
a torque balance (see [21 E] for details). 

In equations ([!])-([3]), we have noted A = A^-t-iA, the spin-over damping factor, with 
A ~ —2.62 -|- 0.259i for the sphere. For P 1, an exact expression of A has actually 
been obtained [H], which undergoes viscous corrections for finite values of P [TH |26] . 

Even if equations ([I])-([3]) are obtained without any inner core, corrections have been 
proposed in the case a = c to take an inner core into account. Using the dimensionless 
inner radius r*, it has been proposed to simply modify A by the factor (1 -|- r/)/{l — rf) 
for a no-slip inner core PI, and by 1/(1 — rf) for a stress-free inner core [36] . 

It is possible to give a geometrical interpretation of equations ([I])-(|3]), which allows 
to obtain analytical constraints on the solutions. Since these constraints are interesting, 
but do not allow to really constrain the range of parameters allowing multiple solutions, 
we give this geometrical interpretation and the associated constraints in Appendix A 


3. Multiple stationary solutions in precessing spheroids 

As shown by Noir and co-workers m, equations 0-0 can actually be recast in a 
unique one, namely equation (24) of [21], which is exactly the implicit solution (3.19) 
of|g. Little algebra shows that this unique equation can be transformed into a lengthy 
polynomial of degree 14 (which reduces to degree 8 in the sphere) for the unknown 17^. 
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(a) (b) 

Figure 2. Evolution of 17^ = (see equation [T]) , with the Poincare number P for 
E = 2.10“® (dashed line) and E = 10“® (solid line), with c/a = 0.97, which gives 
Ar ~ —2.646 and \i ~ 0.306 (inviscid values obtained from the formula of Zhang and 
co-workers [44]). (a) a = 30°. (b) a = 81°. 


This polynomial feature guarantees us to obtain efficiently all the possible solutions. 
However, one shall be careful since all the polynomial roots are not necessarily solutions 
of equations (II])-(l3|). We thus systematically test a posteriori the obtained solutions by 
replacing them into equations ([I|)-([3|). 

3.1. Examples of multiple solutions 

In £gure[2^, solutions of ([I])-([3|) are represented in function of the Poincare number P, 
for a moderate precession angle of a = 30°. This shows that, for large enough Ekman 
numbers (e.g. E = 2.10“® in the figure), equations ([I])-(El) lead to only one solution 
for each value of P. However, when the Ekman number is smaller than a certain 
critical value Emax, certain values of P can lead to multiple solutions (figure [2|i) . In 
this case, we can delineate three branches separated by a cusp point and a point where 
dVt/dP = oo, and we note Pg and Pres the respective associated Poincare numbers (in 
figure Et, Pres < Ps)- Note that, according to Noir and co-workers m, the branch 
between Pg and Preg is unstable, and cannot be physically realized. 

In hgure Eb, we increase the precession angle to a = 81°, leading to a quasi- 
symmetrical problem when P is changed in —P (naturally, the problem is exactly 
symmetrical P O —P for a = 90°). This is reflected by the vertical quasi-symmetry 
in figure El Considering the case E = 10“®, i.e. the case E < Emax, this figure shows 
that four singular points, delineating five branches, can exist when a is larger than 
a certain value aum- Starting from P = —oo, we note the four singular points as 
Pres < Ps < Ps 2 < Pres 2 - Note that this ordering is the opposite when c/a > 1 (prolate 
spheroids), and that Ps 2 , and Pres 2 only exist for large values of a. We thus define 
non-ambiguously Pg as the cusp point existing for any a, Pg 2 as the second cusp point 
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E a 

(a) (b) 

Figure 3. Within the solid lines, i.e. for P G [Pres', ^s] and P G [Ps 2 ', Pres 2 ], equations 
(IT|)-(|31) admit multiple solutions for cja = 0.97, which gives Xr ~ —2.646 and Xi ~ 0.306 
(inviscid values obtained from the formula of Zhang and co-workers |44] 1. Then, the 
two figures show the ranges of P leading to multiple solutions (hatched areas), (a) in 
function of E, for a = 83.7°, and (b) in function of a, for E = 10“®. We clearly see 
that Emax, obtained here around Emax ~ 1.1 • 10“®, is actually given by Pres = Ps, 
whereas aum, obtained here around a^m ~ 79°, is actually given by Pres 2 = Ps 2 
{Pres 2 = Ps 2 also defines Emax 2 , obtained here around Emax 2 « 3.1 • 10“®). 


appearing when a is large enough, and Pres, Pres 2 as the points where dfl/dP = oo. 
The point Pres exist for any a, whereas Pres 2 only exists when a is large enough (in any 
case, PsPres < VsPs < Pl.Ps 2 < P?.Pres 2 )- 

In hgure IH another point of view is proposed on these ranges of P where equations 
©-ii) admit multiple solutions. In hgure [Sti, we hx a = 83.7° > aum, which gives 
two zones with multiple solutions. The zone between Pres and Pg exists as soon as 
E < Emax, with Emax giveu by Ps = Pres- The second zone, between Pres 2 and Ps 2 only 
exists at large precession angle (a > aum), and for E < Emax 2 , where Emax 2 is given by 
Ps 2 = Pres 2 - hi this hgure, we clearly see the seven quantities we are interested in, i.e. 
Pres, Ps, Ps 2 , Pres 2 , Emax, Emax 2 , and aum, which are bouuds of the multiple solutions 
areas. Our goal is thus to obtain estimates of these various quantities, in order to obtain 
analytical insights on these multiple solutions zoness. 

3.2. Busse stationary solutions seen as fixed points of a dynamical system 

First, we did not manage to obtain analytical results directly from equations ©-([3]). 
Second, within the mutiple solutions range of parameters, it would be interesting to test 
if these solutions can be realized experimentally, i.e. to study the stability in time of 
these solutions. To do so, equations ©-© have to be obtained as hxed points of a 
dynamical model for f2. Fortunately, these two issues have been recently tackled |25] . 
Without any approximations, they show that O is governed by (see equations A14-A 
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16 of [25] for a spheroid) 


dQx 

dVty 

dt 

dVtz 


Pz^y ~ e [Pz^y + + PPu ■ ^xi 

Px^z ~ Pz^x + e [Pz^x + + PPu ■ ^yi 


— Px^y — e Px^y + PP V • e 


Z-! 


(4) 

(5) 

( 6 ) 


with the viscous term CTy. Contrary to equations ([I|)-(l3l), this reduced model does not 
assume small flattening for the inviscid part, and the results obtained in the inviscid 
limit E = t) will thus be valid for any oblateness. We will thus pay a particular attention 
to this limit, where this model gives more general and accurate results than equations 

0 - 0 . 


So far, the expression of the viscous term has not been obtained in the general 
case. However, in the limit -C 1, and if the angle between the fluid rotation vector 
and the container rotation vector is small, one can obtain an expression of CTy using 
the linear asymptotic expression of spin-up and of the spin-over mode (see the so-called 
generalized model, given by equation 2.28 of [25]): 






/ 










n 


f Qy \ 


I 'so 


y 


+ x 




sup 




/ \ 

Q,y 

y Ptz j 


,(7) 


with 


A..P = - F (1-1/4,1/2], |3/4], 1 - , (8) 

where F is simply the gamma function and F{n,d,z) is the usual generalized 
hypergeometric function, also known as the Barnes extended hypergeometric function 
(see respectively chap. 6 and 15 of Cl)- 

Naturally, equations da-o, obtained in the very particular limit P 1 and 
e -C 1, are recovered as the hxed points of the dynamical model (|l])-([6]) in this limit. 
For instance, taking (jlJixHj;-!- (I5])xffy-|- (E^xfl^ yields 

fPQQ 

in~k)-n= ^ (9) 


Then, in the limit ePx/y/E 1, we recover the so-called no spin-up condition ([T]). This 
condition is thus not valid in general for a spheroid of arbitrary ellipticity. In this limit, 
which is the limit of validity of equations ([I])-(|3]), the viscous term PFj, reduces to 




f Plx \ 

A* 

+ — 
n 

/ 

fly 



CT^ = {EnflP 


fly 


fix 





- 1 y 

1 

0 

J 



( 10 ) 


by simply using the no spin-up equation 

Finally, ([I])-([3]) are simply the hxed points of the dynamical system given by 
equations (|l])-([6]) and (na, and the solutions stability can be calculated in this 
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framework. However, we still need a simpler set of equations to perform tractable 
analytical calculations. To do so, we follow |25] , who showed that the viscous term flTOl) 
can be very well approximated by the so-called reduced form (equation 2.29 of their 
work [25]) 


CT, 


/ 

x^/e 


12 3, 
Xly 


\ 


V -1 / 


( 11 ) 


which is linear in f2. One can notice that the linear viscous terms of the reduced model 
does not include the coefficient Xi, and thus neglect its influence (typically, a small 
viscous modihcation of the values of P where the Poincare flow undergoes a resonance, 
see [2S] for details). 

The analytical calculations presented in this work have been performed with the 
computer algebra system MAPLE. They have been compared with numerical solutions 
of algebraic and differential equations (e.g. equations [T]l3] and mini respectively) solved 
with the dedicated MATLAB script FLIPPER. This home-made script allows to solve 
efficiently all the equations described above, either by time-stepping or by directly 
looking for all the possible steady solutions. This script can also solve the system 
of equations proposed by [S] and [2S] for precessing triaxial ellipsoids, implementing the 
various viscous terms (equations ([71 [TOl and[TT|), in each case. The script FLIPPER and 
its documentation are provided in supplementary materials. 


3.3. Analytical estimates, using the so-called reduced model 


The calculation details are given in Appendix B, and we thus only report below the 
important analytical results and steps of this calculation. 

Focusing on stationary solutions of equations (Il|)-(l6|) with the viscous term flTT]) . 
these equations can be recast in a unique polynomial of degree 3 for the unknown 
which allows tractable algebra investigations. For the sphere (e = 0), this polynomial 
reduces to a linear polynomial, and the explicit solution for is then 

sin 2a 


2(p2 + A2p)’ 


Q.y — 


Xlz = 


Psin(a) XrVE 
P2 + A2 P ’ 
XlE + P2 cos2 a 


( 12 ) 

(13) 

P2 + A2P • 

For the sphere, there is thus always a unique solution for the reduced model. Note 
also that the solution ffT^ - flTT|) is interesting in the limit P = 0. Indeed, calculating a 
uniform vorticity solution of the Euler equations for a precessing spheroid leads to the 
so-called Poincare flow, which has a free parameter (as a consequence of the inviscid 
hypothesis). Among this class of solutions, the one usually chosen in the literature is 
dehned arbitrarily by putting 122 = 1 [35]. Here, equations flT^ - ffTTll show that, for 
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the sphere, the solution in the limit of vanishing viscosity is actually the Poincare flow 
dehned with = cos^a. 

We consider below the case of a spheroid (e 7 ^ 0). Little algebra shows that the 
multiple solutions boundaries Ps and Ps 2 (see hgure [3]) can be described by a unique 
boundary Eg = f{P), where f{P) is an analytically known root of a polynomial of 
degree 3. An expansion for P 1 gives the following tractable expression at order 
0{P^) 


X^Eg = 


— slcoP^ + , 2 'Si[4Co(l + e^) + (1 + e) — Sec^jP"^, (15) 


e ■ 4e^ 

where the branch Pg is obtained for sign(e)P > 0, and the branch Pg 2 for sign(e)P < 
sign(e)P^* 2 ’^. In the planetary relevant limit e 1, equation flTSl) gives 

A2/3 


P = TTiTvS e''" L 

S,- Co' ^ 


(16) 


At the order 0{P^) of equation flT^ . putting = 0 leads to two solutions for P. The 
hrst solution is P = 0, corresponding to the inviscid limit P'™ = 0 of P^. The second 
solution corresponds to the inviscid limit P*™ of Ps 2 , given by 


ryinv _ 

s2 “ 


4e(e - l)co 


l + e + c2[3-e(9-4e)]' 

Note that a compact accurate expansion of flT^ can be obtained at the order (P(P^°) 
for the particular case a = 7r/2: 


11 


91 


51 


969 


— Cl + C2 + xCs + —rCi + TxCs + -rCe + “XTrCr + 


4807 


Cs, 


(18) 


2’"' 4^ 16^ 4^ 32^ 64 

with Cfc = (1 + Note also that the leading order of equation ([15]) 

vanishes for a = 7 r/ 2 , which imposes to consider the next order. In the limit e 1, 
equation (IT^ then gives 


P = E'/* + 0(e®''^) for o = 4 


(19) 


Since Pres and Pres 2 are actually weakly dependent on E (see e.g. £g. [3]), it turns 
out that their inviscid limit values, noted respectively Pfl^^ and Pfl^g 2 -i provide good 
estimates of these quantities. They are given by 

3e[2(l + 7e[2(l — 


pinv ^ 
^ res ^ 


and 


a 1 : 


X = 7r/2 — a <C 1 : 


pinv ^ 
^res2 ^ 


2(1 - e)5/3 


+ 


4(1-e)3 


+ (P(W) 


e(l — e) e(3e^ —5e + 4) 


( 20 ) 


2 \/l + e 2(1 + e) 


X + 


4(1 + e)3/2 


x^ + 0{x^) 


e(l — e) e(3e^ —5e + 4) 


2 \/l + e 2(1 + e) 


X 


4(1 + e)3/2 


x"^ + 0{x^), 


( 21 ) 


for x = 7 r /2 — a -C 1. It is naturally satisfying that the leading order of Preg for a -C 1 
recovers the usual linear inviscid resonance P^ = e/(l — e) of the Poincare flow [2^ IB]. 
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(a) 



(b) 


Figure 4. We fix here c/a = 0.97. (a) Thick solid lines correspond, from bottom 
to top, to P™//, P/™, P/™, P/es 2 (roots of Dq, See [Appendix B| for details), whereas 
the dashed lines are the associated expansions, respectively given by equation (120]) . 
pmv _ equation ([TTl) . and equation ([2T]). The vertical dashed line is given 
by equation (1221) for (b) Thick solid lines corresponds to Emax and Ejnax 2 , 

respectively estimated by the two intersection points between Pg and as well 

as Ps 2 and (be. the solutions of the system A = 0, Dp = 0, see [Appendix Bj 

for details). The dashed lines are the assiociated expansions, respectively given by 
equations (1^ and (HI- We have used here = —2.62. 


It is important to note that we provide here, for the hrst time, an analytical estimate 
of the higher-order corrections to the linear resonance Pr of the Poincare flow, i.e an 
estimate of the so-called non-linear resonance [23]. It is also satisfying to retrieve the 
quasi-symmetry of the problem (with respect to P = 0) in the expressions of Pres and 
Pres 2 for a: = 7r/2 — a -C 1. 

Based on these estimates, one can calculate the critical values aum and Emax- The 
critical angle aum is well estimated by its inviscid limit given by 


= arcos 


( / ^ + " ^ 
\^V 27e2 - 53e + 28 j ’ 


( 22 ) 


which is decreasing between 7r/2, reached for e = — 1, and 0, reached for e = 1. Finally, 
Emax (resp. Emax 2 ) is simply estimated by the intersection point between P/Z^g and P*"''“ 
(resp. P/eh and P* 2 ^). We thus obtain 




^ ^ J , 3(2^3-miar'” _ (13,/3-18)Ka)" ^ 


X = 7r/2 — a <C 1 : 


2^/3 2 

-8 + 50 ^ (5 - 30 ) a ; (39 - 230 ) a;2 


( 23 ) 


+ 0{x^ 


4 ' 2^ 10^2 

with .^ = (1 -b e)^/^/(l — e), and the golden ratio 0 = (1 -b \fS)/2. Similarly, Emax 2 is 
given by 

^ Emax2 


-8 + 50 (5 - 30 )^ (39 - 230)^2 , 38 


2^ 


-X 


loe 


-x^ + 0{xfl, 


( 24 ) 
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(a) (b) 


Figure 5. Hatched areas bounded by thick solid lines represent the two multiple 
solutions zones of equations o-®, whereas the dashed lines are our analytical 
estimates. Parameters: (a) small precession angle of a = 2.86°, with c/a = 0.98, 
which gives Xr ~ —2.637 and Xi ~ 0.290 (inviscid values obtained from the formula of 
Zhang and co-workers [H]); (b) large precession angle [a = 83.7°), with c/a = 0.97, 
which gives Xr ~ —2.646 and Xi ~ 0.306 (inviscid values obtained from the formula of 
Zhang and co-workers |44 ) ). 


with X = 7r/2 — a -C 1. The problem symmetry for a = 7r/2 is recovered in the 
expressions (|2^ and fl2Tl) since Emax = Emax 2 in this case. 

These various estimates have been represented in hgure (jl]), showing that they 
capture correctly the multiple solutions zones features in the inviscid limit. Now, we 
can wonder how these estimates compare with the exact solutions of equations (fT |) -(j3 |) . 
In hgure m equations (II])-(IH]) are solved for two different precession angles and spheroids, 
and the ranges of parameters where we have multiple solutions are located between the 
thick solid lines. As expected, for a small precession angle (hgure |5^), there is a unique 
multiple solutions zone, whereas two multiple solutions zones exist for a large precession 
angle (hgure |5 ]d, where we have added to hgure [3^ the various analytical estimates 
we have obtained). Figure |5] conhrms that the analytical estimates obtained from the 
reduced model capture quite well these multiple solutions zones. The previously derived 
expressions allow thus to bound quite accurately these zones, especially in the inviscid 
limit. 

4. Solutions stability 

4 . 1 . Estimates of the Jacobian eigenvalues 

Having localized the ranges of parameters where multiple solutions exist for the how in 
a processing spheroid, one can wonder if it is possible to observe experimentally these 
multiple solutions, i.e. what is the stability of these solutions. According to Noir and co¬ 
workers [2l], the branch between Pg and Pres is unstable and cannot thus be physically 
realized. However, very few details are provided, and we propose thus to reinvestigate 













Bistable flows in precessing spheroids 


12 





(c) (d) 

Figure 6. (a) Time-evolution of ftz (solid lines), starting from the three perturbed 
fixed points (dashed lines) of the dynamical model dH)-®, with the viscous term (fTU|) . 
relevant for equations ((I])-®, for a = 81°, P = —0.011, E = 10“®, and c/a = 0.97, 
which gives Xr ~ —2.646 and Xi ~ 0.306 (inviscid values obtained from the formula of 
Zhang and co-workers [44 ] ). The (unstable) intermediate solution can evolve toward 
the other fixed point with a different initial perturbation (dashed line), (b) Perturbing 
the unstable fixed point in all the directions of the space (flj,, fly, il^), we show here 
that the time-evolution of this solution systematically follows the same path in this 
space (circles correspond to the three fixed points), (c) Same as figure a, with only the 
solution starting from the lowest = 0) (stable fixed point), (d) Same as figure a, 
with only the solution starting from the perturbed unstable equilibrium point. Time- 
evolution of riz is plotted in this way to show the clear exponential growth, with a 
measured growth rate of cr ~ 0.0037. 


this issue here, using the dynamical model described in section 13.21 

We consider the linear stability of the equilibrium solution and we investigate 
the fate of the flow r2° -I- e (where |e| = \(€x,ey,€z)\ 1). Inserting this ansatz in the 

dynamical reduced model (equations HHl] and [TT]) leads to 


-XrVE Pz{l - e) - 
— Pz{l — c) + — Xr\fE Px + oPl^x 

0 —Pxifl + e) —Xy^fp 


dt 


e = MyS, (25) 
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Figure 7. Fixed points of (H])-®, with the viscous term (fTUl) . for a = 81°, E = 10“®, 
and c/a = 0.97, which gives ~ —2.646 and \i ~ 0.306 (inviscid values obtained 
from the formula of Zhang and co-workers [H]). The colorbar shows max^ 77.e(/ii), i.e. 
the maximum of the real part of the eigenvalues, indicating unstable solutions when 
positive, (a) a = 30°. (b) a = 81°. It is clear than the fixed points are unstable only 
on the branch linking to Presj and on the branch linking Pg^ to Pres 2 - 


where the eigenvalues pi (with z = 1, 2, 3) of the Jacobian matrix Mr characterize 
the stability of the equilibrium f2°. Naturally, we can obtain similar Jacobian matrix 
Mg and Mb using respectively the viscous term ([7]) of the generalized model, or its 
simplihed form (1T0|) in the validity limit of equations (II])-(l3]). 

Considering the dynamical model (jlj)-(l 6 ]), with the viscous term flTOll . relevant for 
equations ([I])-(l3]), we show in £gure[ 6 K the time evolution of 12 ^, when we start from the 
three perturbed hxed points (dashed lines). The solution with the intermediate initial 
condition flzfl = 0 ) is clearly unstable, and evolves toward one of the two other hxed 
points, depending on the initial perturbation. As a complementary view, we show in 
hgure [ 6 b the time-evolution of this unstable solution in the space (fl^,, 12 ^,, 12 ^), when 
exploring all the possible perturbations. Note that the time-evolution of this solution 
systematically follows the same path in this space. 

Figure | 6 h is a zoom of hgure [ 6 b, showing only the solution starting from the smallest 
flz{t = 0). This solution comes back to the hxed point as a damped harmonic oscillator, 
decaying exponentially towards the equilibrium point at a decay rate of a ~ —8.9 ■ 10“^, 
and oscillating at the frequency 9.6 ■ 10“^. For this solution, the three eigenvalues pi of 
Mb are pi ~ —0.63 ■ 10“^, p 2 = (—0.88 -|- 9.5i) ■ 10“^, and ps = 7 I 2 , complex conjugate 
of p 2 - The real parts of the three eigenvalues are all negative, conhrming the stability of 
the hxed point, and the eigenvalues are in 1 % agreement with the measured decay rate 
and oscillation frequency. Using the eigenvectors, we have naturally chosen here the 
right initial perturbation to test p 2 , but the three eigenvalues can actually be checked 
with the diherential equations solution, by calculating numerically Mr at t = 0 , using 
de/dt. As expected, this approach gives very close eigenvalues. 

Figure [ 6 b is a zoom of hgure [ 6 b, showing only the solution starting from the unstable 
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fixed point. The solution shows a clear exponential departure from the equilibrium point, 
with a growth rate a ~ 3.7 ■ 10“^, and an oscillation frequency of 7.8 ■ 10“^. For this 
solution, the three eigenvalues of Mb are ~ 3.73 ■ 10 p 2 = (—4.6 + 5.7i) ■ 10 and 
h'S = Wi- Since the real part of one eigenvalue is positive, this hxed point is conhrmed to 
be unstable, and the growth rate is in excellent agreement with the theoretical prediction 
(oscillation frequencies are in less good agreement). 

Finally, the eigenvalues for the last solutions are pi ~ —2 ■ 10“^, /i 2 = (—0.28 + 
2.4i) ■ 10“^, and ps = 7 I 2 , in excellent agreement with the eigenvalues obtained from the 
time-evolution of the solutions. 

In hgure [71 we consider the same parameters as hgure El and we vary the Poincare 
number P, plotting maxj77e(/ij) as a colorbar to indicate the stability of each solution. 
It allows to clearly see that the hxed points are unstable only on the branch linking P^ 
to Pres, and on the branch linking Ps 2 to Pres 2 - In the particular case a < aum, we thus 
recover the conclusion of Noir and co-workers [24j . 

One can also investigate analytically, in particular cases, the eigenvalues of the 
Jacobian matrix. Considering for instance the matrix fl25|l . it is clear that this requires 
an analytical estimate of 0° (with 0°, 0° respectively given by equations IB.21 IB.ip . 
The full analytical solutions (given by the roots of equation IB.SP are very lengthy and 
are thus not tractable, and we will thus work in two simplifying limit cases. The hrst 
limit case is a <C 1 , where 


Vtz = 1 + 


(l + e)(p2_e2) 


P^ + 2 


(l-e2)(2p2_e2) 


+ C>(P^ + K^) (26) 


with K = —Ar's/P- The second limit case is a = 7 r/ 2 , where 

e 


P < Pre, = 


2vTTe ■ 


= _ - _ 

" (l-t-e)P2 




2 


P| < Pres & P < Emax ■ O* 


(l + e)P^ / 

e2 



p2 

(l + e)P2 


(27) 


IPI < P 

\ 1 \ ^ T'i 


0 “ = 1 


l + e 




p2 


at the order 0{P‘^ + K^). In this case, we have clearly three possible solutions, the lower 
one fl[, the unstable one O*, and the upper one flfl 

Based on these expansion, one can estimate the eigenvalues in the limit a <C 1, 


/ii = -K, (28) 

P 2 = — E + i{l — e){P — Pr), (29) 

h3 = F 2 (30) 

where Pr is the (linear) Poincare resonance (see equation |20|) . Here, the solution is 
stable, and, if perturbed, the oscillation frequency |Xm/i 2 | = [Xm/isl is proportional to 
P — Pr- Note that the third eigenvalue is systematically the complex conjugate of p 2 , 
and is thus not considered hereinafter. 
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(a) (b) 


Figure 8. Fixed points of dl])-® for a = 90°, E = 10“®, and c/a = 1.1, which gives 
Xr « —2.536 and Xi ~ 0.120 (inviscid values obtained from the formula of Zhang and 
co-workers [H]). (a) With the viscous term m- The colorbar shows max^ 7?.e(/ii), 
and the dashed lines are given by equation (EZD, obtained with the reduced viscous 
term (fTTl) . (b) Dots are calculated with the viscous term (fTUl) , crosses with the reduced 
viscous torque m, and the dashed line corresponds to equation (1551) . 


In the limit a = 7r/2, we have 

hi = -K, 

= -k + ipVTV 


for fl[, and 


h2 ~ 


hi = 


(31) 

(32) 


-K+[{1 + + eflP^yxY^^ 


(33) 


h2= -K- 


[(1 + e)KP^fl/fll - iv^) [(1 + + i\/3) 


+ 


3e3 


(34) 


for hi! 


/i“= -K+{l + e)KPye‘^, 


(35) 

/i“= -K + i[e-P\l + l/e)], (36) 

for 12“. The solution 12* becomes unstable due to the driving term x = [(1 + e)KP‘^Y^^ 
in equation fl33ll . which is thus the control parameter of this instability. Note that the 
real part of the two other eigenvalues of 12* is always negative, denoting a saddle-node 
bifurcation of 12*. 

In hgure [Hfei, we consider the limit case a = 90°, which is a perfectly symmetrical 
conhguration for P and —P. In £gure[8], the growth rate/decay of the dynamical model 
corresponding to the Busse equation are shown, and compared with the results of the 
reduced model. The estimate fl5^ is also shown, capturing correctly the global behaviour 
of the growth rate/decay, with a maximum reached for a Poincare number of 

,2 7^2 \ 3/8 


P = 

mciT 


3e^K^ 


' ^{l + e)K \ 10 


( 37 ) 
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giving an analytical estimate of the most unstable Poincare number. 

One can wonder if we can interpret the previous results in the framework of the 
dynamical systems theory, and, for instance, obtain typical bifurcation diagrams. In 
section 14.21 we focus on the case a = 7r/2, where exact calculations can be performed 
analytically without any expansion or hypothesis, in order to investigate the onset, i.e. 
the zone around Emax- 


4-2. Nature of the bifurcation 


Considering the cubic equation governing the fixed points of given by equation (IB.511 . 
one can obtain easily obtain P = f{e, K, as roots of a quadratic equation (without 
any hypothesis). Equation dP/dQz = 0, which defines Pres = dflz/dP = oo, is, in the 
general case, a polynomial of degree 6. However, it reduces to a cubic polynomial for 
a = 7r/2, given by 

- 2e^nl + -K^ = 0. (38) 


Note that, for iC = 0, we obtain = 1/2 and = 0, which naturally allows to 
recover P = /(e, K, H^) = e/(2\/l + e) = Pres- The discriminant of the polynomial 
fl35]) is such that 


^res _ 


Kfl27K‘^ - e^) 
108e4 


= 0 


K = -XrE = 


33/2- 


(39) 


Then, when X^E > 127^ equation fl38|) have only one real solution, and, at most, three 

solutions otherwise. Since we are considering the resonant Poincare number P, this 
shows that we have obtained the rigorous value X^Emax = ^{27 for a = 7r/2. At this 
particular point K = e/(3^/^), the three roots of equation (l38|l are equal to = 1/3, 

and P is given by 


JD^max 
^ res 


1 /32 

= f{e,K,nz) = — 


2\ 27 yiT^' 


We now focus on the onset of the instability by perturbing E„ 
Ekman number 

^2 


Em.ax 


27X1 


(l-'S). 


(40) 

into the perturbed 

(41) 


which is equivalent to consider the perturbed quantity K = e(l — 5/2) /We then 
have 


JD^max 
^ res 


= f{e,K,nz) = 


36\/rTe 

When 5 < 0, the unique solution of equation 


\/ m -\/6 5 + \/ 2 ( 5 T 2 ]+ 0 (^ 2 ). ( 42 ) 

is thus 


n 


basic 

z 



(43) 


at leading order in 6 . This corresponds to the basic state. Note that, for 5 < 0, we 
naturally do not have dflz/dP = oo anymore ; this basic state simply corresponds to 
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Figure 9. Results for the fixed points of (jj])-® with the reduced viscous term 
(ED, using c/a = 1.1, a = 7r/2, and the perturbed Ekman number (1411) . such that 
K = e(l —d/2)/(3^'^^). (a) Comparison of the exact Pres (solid line) with the expansion 
(l42)l . (b) Comparison of the two exact solutions, 11“*’ and (solid lines), bifurcating 
from the unique solution (existing for <5 < 0) with their respective expansions 

(HU, and (USD. 


the flow obtained when P is given by equation fH2ll and E by equation fHTll . When 
5 > 0, equation fl38ll has the two following solutions , 


/q 

QUP ^ phasic (44) 

9 

e? + (45) 

at the order (9(5®/^). Note that the two solutions, the upper one and the lower one 
are not symmetrical with respect to the basic state In figure |9l we compare 

the expansions (H2D . (jHD, and (H5ll . with their exact counterparts. The agreement is 
very good, and one can notice the clear bifurcation of the solution Viz = 1/3 into the 
two solutions Vt'fl and hi/’ for 5 > 0. 

One can now calculate how the eigenvalues vary with 5 by using the Jacobian matrix 
(EHD- We obtain, at the order 


e 


up 

P2 


up 

P3 


x/3 


e e 

for the solution 0“*’, and 


= —^ + 1 


Pi 


P2 


^/3 


5, 




6 


I 611\/^ + 1656\/3 
69(611 + 72\/^) ^ ) ’ 


\/Q9 


^ 611v^ + 1656\/3 A 
69(611 + 72v^) ^ ) ’ 


(46) 

(47) 


(48) 

(49) 


e 


e 
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Figure 10. (a) Considering the figure 11 of the experiments m in a spheroid for 
a = 7r/2, and a/c = 0.9, which gives \r ~ —2.738 and \i ~ 0.404 (inviscid values 
obtained from [M]), the rounds indicates the two scaling laws, in P ~ and 

P ^ '/E, proposed by [11] for the experimental instability onset (P = 2 for 
the empty rounds, and P = 5.5 ^/E for the solid rounds, determined by fitting their 
experimental results). The solid line represent the multiple solutions zone for equations 
m-m, without any adjustable parameter. According to estimate m, the lower solid 
line scales as when a 7r/2, and as E^!^ when a = 7r/2 (eq. [lEl which do 

not use any approximation, rather give here, probably because Ai ^ 0 in these 
equations), (b) Considering a triaxial ellipsoid (hja = 0.8, cja = 0.7, P = 10“^, 
a = 0.3, P = —0.0255) and using the reduced viscous term (ITTl) . we start from the 
three possible steady solutions for H obtained in the spheroid with the same cja. The 
time-evolutions clearly show that two possible flows can exist, and both are periodic 
in time and remain close from the stable solutions of the spheroid. 


for the solution At the considered order, the solution = 1/3 bifurcates for <5 > 0 
into the marginally stable solution VLfl’ and the stable solution 

To summarize this study at a = 7r/2, the hxed point of the dynamical system (jUlH]) 
with the viscous term flTT]) loses stability when X^Emax < e^/27, as the real eigenvalue 
(H6|l crosses 0, which indicates a pitchfork bifurcation. Besides, one can notice the usual 
square root dependency of the amplitude above the onset (equations HU and HS]) , which 
is clear in hgure [9l The two other eigenvalues are complex congugates with a negative 
real part, indicating that the pitchfork bifurcation originates from a saddle-node. 

5. Discussion 

We can conclude from the previous sections that two stable solutions can coexist on 
the branches linking Pg to Pres, and Ps 2 to Pres 2 - However, to the knowledge of the 
author, these coexistent solutions have never been observed in the literature, neither 
experimentally, nor in numerical simulations. In presence of a strong enough noise, the 
dynamical system could jump intermittently from one stable hxed point to the other. 
In this case, this could have been interpreted as an instability. In hgure [101 'w® present 
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the very recent experimental results obtained by Goto and co-workers m, and we 
compare them to the multiple solutions zones given by equations (II])-Q. A rather good 
agreement is found, without any adjustable parameter. Note that this apparent good 
agreement could also be linked with the fact that, around Pres, fhe solution suddenly 
jumps to another branch. This drastically modihes the flow, and the new flow may then 
excite an instability (e.g. an inertial one, as described by na). 

One can wonder if the Earth has undergone, during its evolution, parameters 
allowing states with possible coexistent flows. Using the values given by EH for a Earth- 
Moon distance varying between its current value and half of this value, and assuming a 
constant flattening equal to its current value e ~ —0.003, we obtain that P remains of 
the order P ^ —10“^, which is larger than Pg G [—10“^; —10“®]. A multiple solutions 
state is thus not expected for the Earth. One can ask the same question for the Moon. 
Based on EH and [7], considering a Earth-Moon distance varying between its current 
value and half of this value, and assuming a constant flattening equal to its current value 
e ~ —2.5 ■ 10“®, we obtain that P remains of the order P ^ —10“^, which is smaller 
than Ri e ~ —2.5 • 10“®. Contrary to the Earth, the Moon has thus a Poincare 
number which is too small for multiple solutions, but one can notice that planetary 
typical values do not allow to discard the possibility of multiple solutions on simple 
orders of magnitude arguments. 

Note that, because of its synchronized state, the Moon is rather a precessing triaxial 
ellipsoid than a spheroid rotating along its symmetry axis |25]. The solutions are then 
time-periodic, and one can thus wonder if such multiple solutions can still exist in this 
case. Figure fTOb shows that the model proposed by [25], and solved by the script 
FLIPPER (supplementary material), allows these multiple solutions when the reduced 
viscous term flTTl) is used. As already noticed by [25], the solutions in the triaxial 
ellipsoid remain close from their analogs in the spheroids. Having checked that these 
multiple solutions also exist using the other viscous terms ([7]) and (ITOll . we thus believe 
that our estimates for the Moon are quite accurate. A more detailed study of the triaxial 
ellipsoid case is beyond the scope of this paper. 

6. Conclusion 

In this work, we investigate the ranges of parameters allowing multiple solutions for 
the flow forced by precession in a spheroid. To do so, we hrst solve the equations very 
efficiently, with the dedicated script FLIPPER, provided as a supplementary material. 
Then, we obtain various analytical results on the solutions. For instance, we obtain 
analytical estimates of the ranges of parameters allowing these multiple solutions, 
and these analytical results are successfully compared with numerical solutions of the 
equations. Finally, the stability of the solutions is analytically obtained, extending the 
results of El in a quantitative manner. This dynamical model approach also allows an 
accurate description of the bifurcation of the unstable flow solution. 

Naturally, it would be interesting to investigate exprimentally or numerically these 
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co-existent solutions, which have not been observed yet. However, the required values of 
the Ekman number are quite small, preventing an easy use of local methods. Moreover, 
usual spherical harmonics codes can only deal with a spherical geometry, where a unique 
solution is always expected. In order to investigate this issue, we plan to develop a 
spectral method designed to deal with spheroids. 
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Appendix A. Geometrical interpretation, constraints on the solutions 

In the set of equations ([I])-®, for a given Hj,, each equation admits solutions in the 
plane {Qx,^y)- In this plane, the solution of the complete system is then given by the 
intersection points of these different solutions locations. 

The so-called no spin-up equation ([1]) describes a sphere of center (0,0,1/2) and 
radius 1/2. For a given solutions are thus on circles of center (0, 0), and radius 

r = 7^,(1-a) = ^1/4 - (H, - 1/2)2. (A.l) 

Since, = 1/4 — — 1/2)^ > 0, we see that G [0; 1], and Hj,, VLy G [0; 1/2]. In 

the following, we are then considering geometrical constraints for a given i.e. in the 
plane 

Equation ([2]) can be rewritten as VLy = fiVLx + 5'i, describing a ligne. This line goes 
through the origin [Py = 0), which leads to two points of intersection with the circle. 
The distance d of this line from the circles centers described by equation ([T]), i.e. to the 
origin, is given hy d"^ = gf/{1 + f^). A solution of the system of equations ([I]12]) is thus 
only possible if an intersection point exists, i.e. if d? < gl < ^^(l — fl 2 )(l -|- //). 
For a; = 1 — 1, the expansion of H^(l — Hz)(l -|- fD/gl leads to < 1, using 

E -C 1. This is not a supplementary constraint on 

Equation (|3]) can be rewritten as VLy = f 2 ^x + g 2 , which describes a line. This line 
is horizontal when B = 0, since we then have Qy = g 2 . This leads to the constraint 
92 — 1 ^ 2(1 “ I^ 2 )(l + /D- For X = 1 — <C 1, the expansion of g 2 /[^z{^ — flz)(l + /!)] 
leads to > 1 — H//(A^E), which is trivial since E -C 1. This is thus always verihed 
in this limit. 

The two lines are crossing when VLy = fiVL^ + gi = f2^x + 92, be. = 

{92 - 9i)/{fi - f 2 )- We thus have |Ha;| < r, i.e. 

(A.2) 

For X = 1 — 112 -C 1, the expansion leads to 

p 2 

_ ■^_X _ 

{Pz + h3 - h2 + Ai\/E)2’ 


> 1 - 


(A.3) 
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Focusing on stationary solutions of equations (|11)-(I6]), these equations can be recast in 
a polynomial of degree 3 for the unknown 12^. Indeed, equation ([6]) gives 


Vty 


( 12 , - l)\r^ 


Pfll + e) ’ 
which can be replaced in equation (|1]), leading to 
[e(12, + P,)-P,)] (12,-1) 

-Px(l + e) 

Finally, using these two expressions, (E]) can be written as 


(B.l) 


(B.2) 


e" 12^ + e(2ec,P - 2c,P - e) 12^ + {P^sl + e^c^P^ + - 2eclP^ (B.3) 

- 2e^CoP + K^ + 2eCoP + clP^) 12, - P' + [e(2 - e) - Ijc'P^ (B.4) 
= 0, (B.5) 


where Co = cos a, Si = sin a and K = Xry/E. One can hrst notice that, for the sphere 
(e = 0), this equation reduces to a linear polynomial, which gives the solution (1T^ - (IT^ . 

We are interested by the number of solutions of equation fIB.Sp . which implies to 
study the sign of the discriminant A of this cubic equation. It turns out that the 
discriminant A is given by an equation of degree 3 in k = A^P, which allows to obtain 
explicitly the roots fci, k 2 , and k^, of A. An expansion for P 1 gives at order 0{P^) 


ki = 


-s^CoP^ + 


1 + e 


and, at order 0 {P^), 

k 2 = - + 2e[co(l - e) - d]P + [d 2 cl 

, „Q1D , LQ „2 


5^(1+ e) -8ec^]P^ 

(B.6) 

e)'^- ^5^(1 + e)]P^ 

(B.7) 

+ e)]P^ 

(B.8) 


where'd = Si^j2{l + e), and ^2 = 2 — — 1. The root ki corresponds actually to P^ 

and Ps 2 (equation [T5|l . and ki is thus important for the caracterization of the multiple 
soution zone. A more accurate expansion of this important root has been obtained for 
a = 7r/2, which is given by equation flT^ . 

We consider now the inviscid limit P = 0, which gives ki = 0. At the order of 
equation flT^ . /ci = 0 for two different values of P. The hrst one is P = 0, which 
corresponds to the inviscid limit P]"''^ = 0 of P^. The second one is given by equation 
031. which is the inviscid limit P^^^ of Ps 2 - 

Focusing on the inviscid limit E = k/X^ = 0, the discriminant A reduces to 
Dq = A{k = 0). Noting that P = 0 is a solution, we consider Dq/P^, which is an 
algebraic equation of degree 3 in P. Since A is the discriminant of equation fIB.Sp . 
the roots of Dq naturally correspond to the multiple solutions zones boundaries in the 
inviscid limit. Thus, the roots of Dq will allow to obtain P,!™, P^eX^ pinv^ 

The solution P = 0 of Pq is naturally the inviscid limit P^*™ = 0 of Ps 2 (which has 
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been seen above). The three remaining roots, which are the analytically known roots of 
Dq/P^, naturally correspond to Pres 2 (equations [ 20 ] and [ 21 ]), and PI™ ^ given by 

the following expansion 


■pinv 

^s2 


4e(e — 1) 
1 + e 


X + 


2e(l-e)(24e2-53e + 19) 3 
3(l + e)2 ^ 


+ 0{x^), 


(B.9) 


where a = tt/ 2 —x <C 1. Note that these expressions of Ps 2 and Pfl^s 2 when a 

is lower than a certain value aum, determined below, and their expansions around a = 0 
are thus not relevant. An expansion of (ITT)) around a: = 0 allows to recover exactly the 
expression flB.9D at order 4 (but the term in x^ differs in the two expansions). In the 
invscid limit, PI™ is thus given accurately by flTT]) for P 1 and arbitrary a, and by 
flB.9p for arbitrary P but 7r/2 — a <C 1. 


The sign of A = 0 gives the number of solution for i.e. the zones where mutiple 
solutions are possible. As shown in hgure [3^, two zones of multiple solutions can exist 
in the plane (P, P). In the inviscid limit, the number of zones is thus directy given by 
the number of solution of Dq = 0 , i.e. by the sign of the discriminant Aq of Dq/P^. 
This discriminant is given by 

4e®(l + ef [(27e2 - 53e + 2fl)cl - 1 - e]^ 


432 [e(e - 3)c2 

which is equal to Aq = 0 for a = 0 or 


1 ]^ 


cos a = ±1 


1 + e 


27e2 - 53e + 


(B.IO) 


(B.ll) 


which naturally gives the inviscid limit a]™) of aum (equation [ 22 ]). 

By dehnition, Emax is the intersection point of Ps and Pres, whereas Emax 2 is the 
intersection point of Ps 2 and Pres 2 - One can thus use our previous estimates to calculate 
them. To do so, we replace P in A = 0 (which dehnes Pg and Ps 2 ), by the inviscid 
estimates P^™ and Pr ™2 of, respectively. Pres and Pres 2 (equations [ 20 ] and [ 2 T]) . This 
gives respectively Emax and Emax 2 (from the roots k = X^E of A). 
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